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SECONDARY  RADIATION-DRIVEN 
THERMAL  INSTABILITY 

I.  INTRODUCTION 

Oran,  Mariska,  and  Boris  (1982)  investigated  the  stability  of  a  one- dimensional 
plasma  at  temperatures  and  densities  typical  of  the  solax  transition  region  and 
corona,  using  both  a  linear  analysis  and  non-linear  time- dependent  numerical  sim¬ 
ulations.  Under  constant  pressure  conditions  with  a  steady  heat  source  H,  some 
regions  of  a  plasma  cool  while  other  regions  are  heated.  The  rate  at  which  this 
process  occurs  depends  on  the  balance  of  H  and  the  radiative  loss  rate  5.  The 
radiation  loss  rate  of  a  plasma,  5,  varies  in  a  complex  way  with  the  mass  density 
and  temperature.  When  5  increases  faster  with  the  mass  density  than  it  decreases 
with  a  lowering  of  the  temperature,  the  cooling  regions  will  tend  to  condense  as 
the  plasma  attempts  to  maintain  pressure  balance.  This  phenomenon  is  termed 
the  condensational  instability.  On  the  other  hand,  the  hotter  regions  will  expand 
and  tend  to  grow  hotter  as  the  steady  heating  (per  unit  volume)  becomes  domi¬ 
nant  over  the  radiation  cooling  process.  This  runaway  temperature  effect  is  usually 
called  the  radiation-driven  thermal  instability  (Priest  19S2).  We  will  use  this  latter 
designation  to  refer  to  either  process. 

Oran,  Mariska  and  Boris  started  from  an  isotropic  static  equilibrium  state. 
They  perturbed  this  state  with  single  mode  disturbances  in  the  velocity  field,  and 
observed  that  the  plasma  condensed  into  one  or  more  cool  dense  regions  with  a  hot 
tenuous  surrounding.  The  condensations  appeared  to  reach  a  secondary  equilibrium. 
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but  they  emphasized  that  they  had  performed  only  one-dimensional  simulations  and 
questioned  whether  their  final  state  might  be  unstable  in  two  or  three  dimensions. 

We  report  here  on  preliminary  results  of  an  investigation  into  the  radiation- 
driven  thermal  instability  of  a  one-dimensional  secondary  equilibrium.  Starting 
from  a  one-dimensional  secondary  equilibrium  state,  we  have  performed  a  series  of 
non-linear  simulations  in  two-dimensions  and  found  that  the  secondary  equilibrium 
is  indeed  unstable  to  random  two-dimensional  perturbations. 

In  §  II  we  describe  the  governing  equations  and  numerical  method.  The  com¬ 
putational  procedure,  which  is  based  on  the  FAST2D  hydrodynamic  algorithm,  is 
the  same  as  that  employed  by  Dahlburg,  et  al.  ( 1987) .  In  §  III  we  describe  the  one¬ 
dimensional  simulation  performed  to  obtain  the  secondary  equilibrium  state.  We 
augment  the  description  of  this  system  given  by  Oran,  Mariska  and  Boris  ( 19S2)  by 
describing  the  energetics  of  the  system.  In  §  IV  the  results  of  the  two-dimensional 
runs  are  described.  We  typically  see  the  system  move  rapidly  away  from  the  one¬ 
dimensional  secondary  equilibrium  state.  In  §  V  we  discuss  rhe  results  in  the  context 
of  the  upper  solar  atmosphere.  In  the  appendix  we  formulate  an  eigenvalue  problem 
for  infinitesimal  disturbances  of  the  secondary  equilibrium. 

II.  FORMULATION  OF  THE  PROBLEM 

We  consider  a  compressible,  two-dimensional,  uniform,  fluid  medium  in  Carte¬ 
sian  geometry  (see  Field  1965,  §  II).  The  evolution  of  this  fluid  is  governed  by 
the  equations  of  continuity,  motion,  energy,  and  state,  which  are  written  below  in 
conservation  form: 


dp 

dt 


+  V  •  (/>v)  =  0. 


i  1  </ 1 


2 


^  +  V-  [(f?  +  P)v]  =V-(KVD  +  ff-5(p,T), 


(1  c) 


P  =  (7  -  1)^,  (1<0 

where:  p(x,y,f)  =  mass  density,  v(x,y,<)  =  (u,  u,  0)  =  flow  velocity,  U(x,y,t) 
=  internal  energy  density,  P(x,  y,  t)  =  total  energy  density  =  U(x,y,t)  +  |/9|v|2, 
P(x,  y,t)  =  mechanical  pressure,  T(x,  y,  t)  =  thermodynamic  temperature,  n(p  ,  T) 
ss  thermal  conductivity  tensor,  H  =  heating  rate,  S(p,  T )  =  radiative  loss  function, 
and  7  =  ratio  of  specific  heats  =  5/3  for  an  ideal  monatomic  gas.  For  the  radiative 
loss  function,  S(p ,  T),  we  use  the  model  employed  by  Oran,  Mariska,  and  Boris 
(1982).  Periodic  boundary  conditions  are  imposed  in  both  spatial  directions. 

The  primary  equilibrium  state  for  the  system  of  equations  la  -  Id  is  given  by: 
p  =  p0  =  constant,  v  =  0,  U  =  t'o  =  constant,  and  H  =  S{pQ,  T0)  =  constant. 
In  this  report,  we  use  the  same  values  for  these  constants  as  Dahlburg,  zt  ai.  viz.. 

( 1986)  po  =  3.07  x  10~14  g  cm-1,  =  4.S5  erg  cm'1,  Hu  =  0.0272  erg  cm-1  s-1. 

k  =  1.14  x  10-6  T2  5  erg  cm-1  s_I  K_l,  and  T0  =  7.29  x  10s  K.  We  use  this  static 
equilibrium  state  to  initialize  our  one-dimensional  code. We  then  perturb  this  state 
and  evolve  the  system  in  time  until  a  state  of  secondary  equilibrium  is  achieved. 

The  numerical  method  used  to  solve  equations  la  -  Id  is  described  in  detail  by 
Dahlburg  et  al.  (1987) .  Here  we  give  a  brief  description  of  the  algorithm.  Centered 
differences  are  used  to  discretize  both  spatial  directions.  For  the  two-dimensional 
simulations  reported  in  this  paper,  the  system  is  a  square  region  9.19  x  10s  cm 
on  edge.  This  size  represents  one  wavelength  of  the  characteristic  perturbation 
used  by  Oran.  Mariska.  and  Boris  (19S2)  for  the  temperature  and  mass  conditions 
under  study.  The  region  is  divided  into  a  100  by  100  grid  with  uniformly  spaced 
intervals  of  9.19  x  10rt  cm.  A  time-step  splitting  scheme  is  used  which  separates 
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the  hydrodynamic  processes  from  the  thermodynamic  processes.  The  hydrodynamic 
terms  are  advanced  in  time  by  the  flux-corrected  transport  method,  with  a  predictor- 
corrector  method  employed  to  time-center  the  source  terms.  For  time  advancing  the 
thermodynamic  terms,  an  iterative  alternating  direction  implicit  scheme  is  used. 
III.  ONE-DIMENSIONAL  RESULTS 

In  this  section  we  describe  the  results  of  the  one-dimensional  run  performed 
to  obtain  the  secondary  equilibrium  state.  To  run  the  two-dimensional  code  in  this 
one- dimensional  mode  we  use  4  uniform  cells  in  the  y  direction,  with  Ay  =  Ax. 
No  variation  in  y  is  permitted.  We  perturb  the  static  equilibrium  state  described 
in  §  II  with  random  perturbations  in  the  x-velocity  component  and  allow  the  sys¬ 
tem  to  evolve  in  time,  where  the  system  length  is  given  by  Lx  =  9.19  x  108  cm. 
We  terminate  the  run  when  we  are  satisfied  that  a  state  of  secondary  equilibrium 
has  been  achieved,  i.e.,  when  the  system  has  stabilized  at  a  different  state.  Our 
simulation  differs  from  that  of  Oran,  Mariska,  and  Boris  in  two  respects.  First, 
we  use  periodic  boundary  conditions.  Second,  we  use  random  initial  conditions  in 
the  velocity  field  so  that  the  system  can  traverse  as  wide  a  range  of  phase  space 
as  possible.  Although  our  simulation  differs  in  some  respects  from  those  of  Oran. 
Mariska,  and  Boris  (1982),  our  results  appear  to  confirm  theirs. 

The  evolution  of  the  system  as  a  whole  is  represented  by  the  internal  energy. 

Figure  1  shows  the  internal  energy,  J  U  dx,  as  a  function  of  time.  This,  and  all  sub- 

o 

sequent  global  quantities,  is  defined  per  unit  length  in  the  r  direction.  The  internal 
energy  remains  approximately  constant,  near  the  primary  equilibrium  value,  as  long 
as  the  perturbations  remain  small.  When  the  perturbations  attain  finite  amplitude, 
as  evidenced  by  the  rapid  growth  in  the  kinetic  energy,  the  internal  energy  drops 
precipitously.  Some  growth  due  to  reheating  is  seen,  and  then  the  internal  energy 
remains  approximately  constant  for  approximately  900  seconds.  We  regard  this 
constant  behaviour  of  the  global  quantities  as  sufficient  evidence  that  the  system 


has  attained  a  state  of  secondary  equilibrium.  Further  evidence  that  a  secondary 
equilibrium  state  has  been  achieved  is  shown  in  figure  2,  which  shows  the  energy 

Ijw 

deficit,  Ed  —  j  /[S(p,  T)  —  H]  dx  dt,  as  a  function  of  time.  The  energy  deficit  re- 
o  o 

mains  approximately  constant  for  almost  700  s,  showing  that  the  heat  energy  input 
and  the  radiation  enrgy  losses  are  about  equal  in  the  vicinity  of  the  primary  equi¬ 
librium.  This  equilibrium  is  destroyed  by  the  formation  of  condensates,  which  are 
extremely  effective  radiators  during  their  formation  phase.  After  some  time  a  state 
of  secondary  equilibrium  is  achieved  in  which  the  spatially  varying  radiation  energy 
losses  equal  the  heat  energy  input.  The  maintenance  of  this  equilibrium  depends 
on  thermal  conduction  and  convection  for  the  transport  of  heat  energy  to  those 
locations  which  are  radiating  energy  most  efficiently. 

Figure  3  shows  the  mass  density  profile,  p(x)  at  secondary  equilibrium.  Figure 
•4  shows  the  corresponding  temperature  profile  T(x).  We  see  that  after  saturation, 
a  cool,  dense  region  is  able  to  exist  in  equilibrium  with  a  hot,  rare  region.  The  max¬ 
imum  mass  density  in  the  condensed  region  is  4.5  x  10~12  g  cm-3,  at  a  temperature 
of  4.2  x  104  Iv.  The  minimum  mass  density  in  the  rarefied  region  is  1.1  x  10~13  g 
cm-3,  at  a  temperature  of  1.7  x  106  K. 

IV.  TWO-DIMENSIONAL  RESULTS 

In  this  section  we  determine  numerically  whether  the  one-dimensional  equilib¬ 
rium  state  described  in  §  III  is  subject  in  two-dimensions  to  secondary  radiation- 
driven  thermal  instabilities.  We  choose  LXJ  =  Lz. 

We  here  describe  one  run  which  is  representative  of  those  which  we  performed 
in  detail.  Figure  5  is  a  vector  plot  of  the  initial  velocity  field  perturbation.  To 
avoid  cluttering  the  plot,  the  velocity  field  is  represented  on  a  reduced  set  of  points. 
The  perturbed  velocity  field  is  determined  by  the  method  used  by  Dahlburg.  e.t 

Lx  Ly  ^ 

al.  (1987).  Figure  6  shows  the  kinetic  energy,  j  f  f  p  |  v  ]2  dx  dy,  as  a  function 

‘oo 

of  time.  Starting  from  an  initial  value  of  about  5.2  x  10lh  ergs  there  is  at  first  a 


very  rapid  drop.  Then  after  a  slight  rise,  the  kinetic  energy  drops  for  a  period  of 
about  500  s  as  it  flows  into  other  forms  of  energy  available  to  the  system.  The  rate 
of  decay  decreases  steadily  and  the  kinetic  energy  appears  to  approach  a  limiting 
value  of  about  1.1  x  1016  ergs. 

L,  L, 

Figure  7  shows  the  interned  energy,  /  f  U  dx  dy ,  as  a  function  of  time.  Starting 

o  o 

from  as  initial  value  of  about  3.45  x  1018  ergs,  this  is  seen  to  decrease  at  an  almost 

constant  rate  for  about  180  s.  Subsequently,  the  rate  of  decay  decreases.  At  the 

end  of  the  simulation,  t  =  727  s,  the  internal  energy  is  equal  to  2.2  x  1018  ergs. 

L*  tm 

Figure  8  shows  the  energy  deficit,  Ed  =  f  f  f  [5(p,  T)  —  H\dx  dy  dt,  as  a  func- 

ooo 

tion  of  time.  This  represents  the  radiation  energy  reduced  by  the  amount  of  heating 
energy.  Initially  this  quantity  is  zero,  representing  the  thermodynamic  equilibrium 
of  radiation  energy  losses  and  constant  heat  energy  input.  It  rises  steadily  from 
zero  at  a  rapid  rate  for  about  180  s.  After  this  time,  the  rate  of  increase  appears 
to  decline.  At  the  final  time  the  energy  deficit  is  equal  to  1.3  x  1018  ergs.  The 
rate  of  increase  of  the  energy  deficit  is  consistent  with  the  rate  of  decrease  of  the 
internal  energy.  This  suggests  that  radiation  is  produced  mainly  at  the  expense  of 
the  internal  energy,  and  that  the  kinetic  energy  is  maintained  by  the  heating  source. 
This  also  suggests  that  the  most  prominent  effect  produced  by  the  perturbing  of 
the  system  is  the  production  of  more  radiated  energy. 

Figure  9  is  a  contour  plot  of  the  mass  density  at  the  final  time.  It  shows 
considerable  concentration  of  mass  in  2  regions  of  the  plot.  The  mass  density 
appears  to  be  kinking,  with  high  density  regions  forming  at  the  extremities  of  the 
kink.  The  maximum  mass  density  is  about  8.6  x  10-13  g  cm-3,  while  the  minimum 
value  is  6.2  x  10-15  g  cm-3.  These  values  differ  by  a  factor  of  about  140.  Figure  10 
is  a  contour  plot  of  the  temperature  at  the  same  time.  The  kinking  is  also  apparent 
in  the  temperature  field,  but  it  appears  to  retain  the  bifurcated  character  of  the 
secondary  equilibrium.  The  maximum  temperature  is  about  1.6  x  106  Iv,  while  the 
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minimum  value  is  1.2  x  104  K.  These  values  differ  by  about  the  same  amount  as  the 
mass  density  maximum  and  minimum.  Figure  11  is  a  vector  plot  of  the  final  velocity 
field.  A  very  active  flow  region  is  visible  at  the  lower  right.  The  twisting  action 
implies  high  vorticity,  which  should  promote  the  formation  of  a  two-dimensional 
state.  Figure  12  is  a  contour  plot  of  the  dilation  factor,  or  divergence,  of  the  fluid. 
A  positive  dilation  factor  signifies  fluid  rarefaction,  and  a  negative  dilation  factor 
signifies  fluid  compression.  This  figure  shows  that  a  rarefaction  zone  develops  ahead 
of  the  high  mass  density  regions,  while  a  condensing  region  forms  behind  this  zone. 

When  very  large  amplitude  perturbations  were  attempted,  say  with  25  times 
the  energy  of  those  described,  the  time-step,  which  is  adjusted  for  each  step  of  the 
simulation  in  response  to  the  hydrodynamic  and  thermal  conduction  time  scales  ( see 
Dahlburg  et  al.  1987) ,  became  so  small  that  even  after  5000  steps  of  the  simulation 
the  total  elapsed  time  was  less  than  half  a  second.  Presumably  in  this  case  the  high 
characteristic  velocity  produced  results  in  an  extremely  short  hydrodynamic  time 
scale. 

V.  DISCUSSION 

In  this  investigation  we  have  seen  the  radiation-driven  thermally  unstable  fluid 
exhibit  a  tendency  to  move  away  from  a  one-dimensional  state  toward  a  two- 
dimensional  state.  In  the  absence  of  any  restoring  force,  it  is  clear  that  the  fluid 
will  not  return  to  its  initial  state.  Unfortunately,  due  to  the  high  cost  of  these  nu¬ 
merical  simulations  we  were  not  able  to  extend  our  calculations  farther  in  time  and 
so  determine  the  long-time  state  of  the  fluid.  The  high  cost  results  from  the  severe 
time-step  restriction  imposed  by  the  thermal  conduction  time-scale.  We  conjecture 
that  eventually  the  sheet  will  break  up. 

Analysis  of  this  phenomenon  is  difficult  due  to  the  high  order  and  extreme  non¬ 
linearity  of  the  governing  system  of  partial  differential  equations.  In  the  appendix 


we  take  a  first  step  in  analysis  by  formulating  an  eigenvalue  problem  for  the  linear 
stability  of  the  secondary  equilibrium. 

In  practice  it  is  difficult  to  observe  the  evolution  of  un.  table  fluids  in  the  solar 
corona  and  transition  region  because  it  has  not  been  possible  yet  to  achieve  the 
necessary  spatial  and  temporal  resolution  in  the  observing  instruments.  Recently, 
faint  superspicules  were  reported  by  Breuckner  et  al.  (1986).  These  superspicules 
were  seen  on  the  disk  of  the  sun  close  to  the  limb,  and  they  appeared  to  form  in  space 
rather  than  to  have  been  ejected  upward  from  the  surface.  Their  available  temporal 
resolution,  at  best  20  s,  but  generally  about  60  s,  would  not  permit  Breuckner  et 
al.  to  observe  the  evolution  of  this  short-lived  phenomenon  (as  little  as  180  s). 
They  did  see  the  superspicules  flex  and  bend  in  the  course  of  their  development. 
These  superspicules  could  conceivably  condense  out  of  the  corona  as  a  result  of 
some  velocity  perturbation  arising  from  the  turbulent  nature  of  the  surrounding 
medium.  Since  the  corona  is  very  tenuous,  the  amount  of  mass  and  consequently 
the  intensity  of  the  radiation  from  these  condensations  would  be  quite  small.  This 
would  explain  their  very  faint  nature. 

Zirin  (1966)  has  stated  that  “loops  and  coronal  rain  typically  occur  as  the 
aftermath  of  flares  and  surges"’.  It  is  not  inconceivable  that  these  phenomena  also 
could  arise  from  the  velocity  perturbations  produced  by  flares  which  then  produce 
condensations  observed  as  loops  and  coronal  rain. 

The  work  performed  here  can  be  extended  by  changing  the  initial  conditions,  or 
the  magnitude  of  the  perturbations,  or  simply  by  adding  the  force  of  gravity  into  the 
calculation  to  see  the  effect  on  the  results.  A  more  difficult  extension  would  be  to 
include  the  effects  of  magnetic  fields  to  simulate  more  realistically  the  magnetofluids 
in  the  transition  region  and  corona  of  the  sun.  Potentially,  this  program  might  Ik* 
applied  to  a  variety  of  phenomena  in  the  sun  or  elsewhere  in  the  universe  to  aid  in 
the  interpretation  of  the  various  phenomena  which  are  observed  (  see  Field  1965). 
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APPENDIX 


In  this  appendix  we  derive  a  generalized  eigenvalue  problem  which  determines 
the  linear  stability  properties  of  the  nonuniform,  secondary  equilibrium.  We  first 
rewrite  the  governing  equations  la  —  lc,  with  the  pressure  eliminated  by  means  of 
the  equation  of  state  Id: 

^  +  V-tf  =  0,  (Ala) 


^  +  V-(*v)  =  -(7-l  )VU, 


(Alb) 


+  V  ■  (Uv)  =  k0V  •  (T$  VT)  +H- p£(p,T) .  (-41c) 

In  this  formulation  we  time  advance  the  internal  energy,  rather  than  the  total  energy. 
The  variable  vIr  is  the  momentum,  pv,  and  £(p,T)  =  -  PpT)-  ■ 

We  first  linearize  equations  .41a  —  .41c  about  the  secondary  equilibrium,  i.c.. 
we  allow  the  field  variables  to  vary  as: 


P  =  Po(y)  +  Pi(x,y,t), 


(-42a) 


V  =  V!(l,(/,t), 


(-426) 


*  =  *!(X,  t/,<), 


(.42c) 
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V  -  Uo(y)  +  U\(x,  y,  t)  , 


(A2d) 


T  =  To(y)  4-  Ti(x,  y,t) ,  (A2e) 

where  the  subscript  0  indicates  the  secondary  equilibrium  field  values,  and  the 
subscript  1  indicates  the  perturbed  field  values. 

Substituting  equations  A2  into  equations  A1  and  ignoring  terms  which  are 
quadratic  in  the  perturbed  quantities,  we  have  to  first  order: 

^-  +  V-*i=0,  (AZa) 

^  =  -(7-1  Nlq,  U36) 


dux 


dt 


+  V 


tfov, 


(.43c) 


.,V  -  (T0?  VT.)  -  -  *(§).(&*.  -  £») 


do 


The  quantities  and  are  evaluated  for  the  equilibrium  state.  Tin 


T  \  '  p 

perturbed  temperature,  Ti,  has  been  eliminated  by  the  first  order  equation  of  state 


in  the  radiation  loss  term.  We  have  not  evaluated  T\  in  the  thermal  conduction 
term  for  notational  simplification. 


After  some  more  algebra: 


dp 

dt 


d*x  O' P, 


dx 


dy 


( .44  c/ ' 


]  1 


mu 

dt 


(.446) 


^  =  -(7-DiS 


(.446) 


dU  TT  (du  dv\  dUo 
dt  °\dx'Jrdy'  dyV 

(.44c/) 

.»v. 

where  we  have  dropped  the  1  subscripts.  This  is  a  set  of  linear  partial  differential 
equations  for  infinitesimal  disturbances  to  the  secondary  equilibrium. 

We  now  perform  a  normal  mode  analysis,  assuming  that  the  variation  of  the 
perturbed  quantities  is  given  by: 


f(x,y,t)  =  f'(y)e‘ax-“-t .  (.43) 

where  a  is  defined  to  be  the  x  perturbation  wavenumber,  and  u:  is  the  complex 
growth-rate  of  the  disturbance. 

If  we  define  D  =  then  upon  substitution  (and  ignoring  the  primed  super¬ 
scripts)  we  have  the  following  set  of  ordinary  differential  equations: 

-iup  =  -iaVt  -  D'pj,  ,  (,46a) 


-iu'Vx  =  -ia( 7  -1)1'. 


(.4Gb) 
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—  iui'l'y  =  -(7  —  l)DU  , 


—iuU  —  —  Uo(iau  +  Dv) 


(-46c) 


•+  k0(DT$)DT  +  k0TJ(D2T  -  a2T) 


(A6d) 


(d£\  /d£\  (To  T0  \ 

Po\dpJ  tP  p0\dTJ  <>\UoU'  p0Pl) 


,dp/T  \oi/p\u  0  po 

We  have  used  the  the  variables  and  T  to  simplify  the  mathematics.  We 

now  eliminate  them  to  reduce  the  number  of  variables  to  the  number  of  equations. 


iapou  +  PqDv  =  iup  , 


(--17a] 


ia(  7  -  1  )U  =  i^'pou  . 


A7h) 


(7  —  1  )Dly  —  iU.'pol'  . 


.47c) 


[p°(^)r~p0(^)p(^)  +K,J0!/>  +  mr°"  +  [CoD~  iDC°). 


+  h(§)p(^)'fK°f5 


U  =  i„'U 


(.4  7,1) 


In  equation  .47<i  we  use  the  following  operators: 


e  =  a2  TJ  Ci  -(DT^K.D  -(DTj)(DCi 


—  Ttj-  I  (  D~  )  4-  2 1  Di,i  )D  +  S|D’ 


1  .4$a  I 
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-<*2To!C2  +  (DT*)<;2D  +  (Dt}){DC2)  +  r0*  [(D2C2)  +  2(D(2)D  +  C 2D2] , 


(.486) 


where  Ci  =  §  and  (2  =  ft- 

Equations  A7  have  the  form  of  a  generalized  eigenvalue  problem: 


Ax  =  ujBx  . 


The  eigenvector  x  =  (  p  u  v  U  )T ,  where  the  superscript  T  denotes  the  transpose 
of  the  vector  x.  The  eigenvalue  is  u>.  the  complex  growth  rate. 
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Figure  3.  Mass  density  at  secondary  equilibr' 
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Figure  6.  Plot  of  kinetic  energy  vs  time  for  two-dimensional  simulation. 
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Figure  7.  Plot  of  internal  energy  vs 
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Figure  9.  Contour  plot 
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